
/*><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
																				
	Replication code for:	
	William A. Sundstrom, Shelby McIntyre, Gregory A. Baker, and Brian Avants,
	Bearers of Bad News: Heterogeneous Effects of Alternative Front-of-Package 
		Labeling Schemes for Nutritional Information
	Journal of Agribusiness
		
	Stata code by William A. Sundstrom, 7/19/2020
	wsundstrom@scu.edu 

*<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

	This replicates output for Tables 2-6 and Figures 4-7 in the paper
	(Table 2 just for the experimental sample stats), as well as some 
	alternative specifications of the regressions mentioned in the paper.
	See paper for details of data, specifications, and calculations.
	
	To get started:
	Set the folder path in line 40 for your machine, copy this file and data 
	set FL_data.csv into it, and create separate folder for output within it.

*/

*===============================================================================
*    Globals etc
*===============================================================================	

	version 16
	clear all
	cap log close
	set more off
	set matsize 800
	
	* packages to install if necessary
	* ssc install xsvmat
	* ssc install outreg2
		
	* Directory globals
	global foodlabel	"/Users/wsundstrom/Desktop/Dropbox/Research/FOPbackup/FOP/Replication" 
	global output		"$foodlabel/output"
	cd $output

	
*===============================================================================
*    Data 
*===============================================================================	

*** Import

	insheet using "$foodlabel/FL_data.csv", c clear
	dropmiss, force
	destring, replace
	sum
	
*** stack the dinners (reshape)

	g str dinner = "chicken"
	g hr = chicken
	g hrc = chickenc
	tempfile bbb
	save `bbb', replace
	tempfile ccc
	save `ccc', replace
	
	foreach x in lasagna meatloaf stirfry { 
		use `bbb', clear
		replace dinner = "`x'"
		replace hr = `x'
		replace hrc = `x'c
		append using `ccc'
		save `ccc', replace
		}
	
	foreach x in chicken lasagna meatloaf stirfry { 
		drop `x' `x'c
		}
	
	replace dinner = "CH" if dinner=="chicken"
	replace dinner = "LA" if dinner=="lasagna"
	replace dinner = "ML" if dinner=="meatloaf"
	replace dinner = "ST" if dinner=="stirfry"
	tab dinner
	
	save `ccc', replace
	
*** new variables

	* treatment dummies
	g hl = condition=="HL"
	g tl = condition=="TL"
	g no = condition=="NO"

	* treatment factor variable
	g treat = .
	replace treat = 0 if condition=="NO"
	replace treat = 1 if condition=="HL"
	replace treat = 2 if condition=="TL"
	label define treats 0 "NO" 1 "HL" 2 "TL"
	label values treat treats
			
	* change in rating pre to post
	g dhr = hr - hrc
	g dhr_log = ln(hr) - ln(hrc)
	g hr_log = ln(hr)
	g hl_hrc = hl*hrc
	g tl_hrc = tl*hrc
	
	* ordinal versions of dhr
	g dhr_ord3 = dhr
	replace dhr_ord3 = -1 if dhr<=-1
	replace dhr_ord3 = 1 if dhr>=1
	g dhr_ord7 = dhr
	replace dhr_ord7 = -3 if dhr<=-3
	replace dhr_ord7 = 3 if dhr>=3
	
	* absolute changes and directions
	g adhr = abs(dhr)
	g dhr_pos = dhr>0
	g dhr_neg = dhr<0
	g dhr_zero = dhr==0
	* dhr_dir=-1 if negative, 0 if 0, +1 if positive
	g dhr_dir = -1*dhr_neg + dhr_pos
	
	* how far is hrc from the mean for that dinner?
	bysort dinner: egen hrc_m = mean(hrc)
	g hrc_dev = hrc - hrc_m
	
	* dinner dummies
	g chick = dinner=="CH"
	g stir = dinner=="ST"
	g meat = dinner=="ML"
	g lasag = dinner=="LA"	
	
	* shorter names
	g ical = impcalories
	g ifat = impfat 
	g isod = impsodium 
	g isug = impsugar	
	g imptotal = ical+ ifat+ isod+ isug
	egen imptot_med = median(imptotal)
	g imptot_high = imptotal>=imptot_med
	
	* relative and absolute nutritional importance ratings
	foreach x in ical ifat isod isug { 
		g `x'_rel = `x'/imptotal
		g `x'_abs = `x'
		egen `x'_med = median(`x')
		g `x'_high = `x'>=`x'_med
		}
	drop *_med

		
*** weighted nutrition (gamma) several flavors

	* we'll use both absolute and relative importance weights,
	* and actual numbers and L-M-H ratings
	
	/*
			Size	Cals	Fat		Sugar	Sodium
		
		ML	14		680		37		29		990
					M		H		H		H
		CH	8		180		3.5		5		490
					L		L		L		L
		LA	8		400		14		5		670
					L		L		L		M
		ST	12		390		4		18		1190
					L		L		H		H
	*/

	* these are the single gammas: weighted index of overall nutritional "badness"
		
	foreach i in abs rel {
	
		* by construction the numerical version of n (content) will average 1
		
		g g_`i'_num = 4*(ical_`i'*680/(680+180+400+390)+ ifat_`i'*37/(37+3.5+14+4) ///
			+ isod_`i'*990/(990+490+670+1190)+ isug_`i'*29/(29+5+5+18)) if dinner=="ML"
		replace g_`i'_num = 4*(ical_`i'*180/(680+180+400+390)+ ifat_`i'*3.5/(37+3.5+14+4) ///
			+ isod_`i'*490/(990+490+670+1190)+ isug_`i'*5/(29+5+5+18)) if dinner=="CH"
		replace g_`i'_num = 4*(ical_`i'*400/(680+180+400+390)+ ifat_`i'*14/(37+3.5+14+4) ///
			+ isod_`i'*670/(990+490+670+1190)+ isug_`i'*5/(29+5+5+18)) if dinner=="LA"
		replace g_`i'_num = 4*(ical_`i'*390/(680+180+400+390)+ ifat_`i'*4/(37+3.5+14+4) ///
			+ isod_`i'*1190/(990+490+670+1190)+ isug_`i'*18/(29+5+5+18)) if dinner=="ST"
		
		g g_`i'_lmh = ical_`i'*2+ ifat_`i'*3 ///
			+ isod_`i'*3+ isug_`i'*3 if dinner=="ML"
		replace g_`i'_lmh = ical_`i'*1+ ifat_`i'*1 ///
			+ isod_`i'*1+ isug_`i'*1 if dinner=="CH"
		replace g_`i'_lmh = ical_`i'*1+ ifat_`i'*1 ///
			+ isod_`i'*1+ isug_`i'*2 if dinner=="LA"
		replace g_`i'_lmh = ical_`i'*1+ ifat_`i'*1 ///
			+ isod_`i'*3+ isug_`i'*3 if dinner=="ST"
			
		}

	* nutrient content for each dinner: n 
		
	g cal_num = 4*680/(680+180+400+390) if dinner=="ML"
	g fat_num = 4*37/(37+3.5+14+4) if dinner=="ML"
	g sod_num = 4*990/(990+490+670+1190) if dinner=="ML"
	g sug_num = 4*29/(29+5+5+18) if dinner=="ML"
	
	replace cal_num = 4*180/(680+180+400+390) if dinner=="CH"
	replace fat_num = 4*3.5/(37+3.5+14+4) if dinner=="CH"
	replace sod_num = 4*490/(990+490+670+1190) if dinner=="CH"
	replace sug_num = 4*5/(29+5+5+18) if dinner=="CH"
	
	replace cal_num = 4*400/(680+180+400+390) if dinner=="LA"
	replace fat_num = 4*14/(37+3.5+14+4) if dinner=="LA"
	replace sod_num = 4*670/(990+490+670+1190) if dinner=="LA"
	replace sug_num = 4*5/(29+5+5+18) if dinner=="LA"
	
	replace cal_num = 4*390/(680+180+400+390) if dinner=="ST"
	replace fat_num = 4*4/(37+3.5+14+4) if dinner=="ST"
	replace sod_num = 4*1190/(990+490+670+1190) if dinner=="ST"
	replace sug_num = 4*18/(29+5+5+18) if dinner=="ST"

	g cal_lmh = 2 if dinner=="ML"
	replace cal_lmh = 1 if (dinner=="CH" | dinner=="LA" | dinner=="ST")
	g fat_lmh = 3 if dinner=="ML"
	replace fat_lmh = 1 if (dinner=="CH" | dinner=="LA" | dinner=="ST")
	g sod_lmh = 3 if (dinner=="ML" | dinner=="ST")
	replace sod_lmh = 1 if (dinner=="CH" | dinner=="LA")
	g sug_lmh = 3 if (dinner=="ML" | dinner=="ST")
	replace sod_lmh = 1 if (dinner=="CH")
	replace sod_lmh = 2 if (dinner=="LA")
	
	* "gamma" separately by nutrient
	
	foreach i in abs rel {
	foreach j in num lmh {
	foreach k in cal fat sod sug { 
		g g_`i'_`j'_`k' = i`k'_`i' * `k'_`j'	
		}
		}
		}		
						
	* interactions of gammas with treatments
	
	foreach i in abs rel {
	foreach j in num lmh { 
	foreach k in no hl tl { 
		g g_`i'_`j'_`k' = g_`i'_`j' * `k'
		}
		}
		}
		
	foreach i in abs rel {
	foreach j in num lmh { 
	foreach k in cal fat sod sug { 
	foreach m in no hl tl { 
		g g_`i'_`j'_`k'_`m' = g_`i'_`j'_`k' * `m'
		}
		} 
		}
		}
					
*** various controls and covariates

	g female = gender=="Female"
	replace female = . if gender==""
	
	* education dummies from panel, less than high school excluded
	g edprov = panel_edcation
	g hsgrad = edprov=="High school degree"
	g somecoll = edprov=="Some college" 
	g assocbach = edprov=="Associate or bachelor degree" 
	g graddegree = edprov=="Graduate degree" 
	foreach x in hsgrad somecoll assocbach graddegree { 
		replace `x' = . if edprov=="" 
		}
	g colldegree = (assocbach==1 | graddegree==1)
		
	* 45plus seems to be what matters
	g age45plus = (age=="45  to 54" | age=="55 or More") 
	g age55plus = (age=="55 or More") 
	foreach x in age45plus age55plus { 
		replace `x' = . if age=="" 
		} 
	
	* misc
	g south = (region=="South Atlantic" | region=="East South Central" | region=="West South Central") 
		replace south = . if region=="" 
	g nonwhite = white~="White"

	* locals for variable sets	
	local xvars1 "female age45plus colldegree"
	local xvars1a "female x25to34 x35to44 age45to54 age55plus"
	local xdinner "lasag meat stir"
	local xvars2 `xvars1' `xdinner'
	local xvars3 `xvars2' "hsgrad somecoll assocbach graddegree inc25to50 inc50to75 incunder25"
	
	* interactions of xvars2 with treatment 
	foreach y in `xvars2' { 
	foreach k in no hl tl { 
		g `y'_`k' = `k' * `y'
		}
		}
		
	* interactions of xvars1 with hrc and treatment 
	foreach y in `xvars1' { 
	foreach k in hrc hl_hrc tl_hrc { 
		g `y'_`k' = `k' * `y'
		}
		}
		
	* interact importance with hrc and treatment 
	foreach i in abs rel { 
	foreach x in ical ifat isod isug { 
	foreach k in hrc hl_hrc tl_hrc { 
		g `x'_`i'_`k' = `x'_`i' * `k'
		}
		}
		}
							
	* interact importance with treatment and xvars2
	foreach i in abs rel { 
	foreach x in ical ifat isod isug { 
	foreach k in `xvars2' no hl tl { 
		g `x'_`i'_`k' = `x'_`i' * `k'
		}
		}
		}
	
	* interactions of xvars2 and gamma by treatment 
	foreach i in abs rel { 
	foreach j in num lmh { 
	foreach y in `xvars2' {
	foreach k in no hl tl { 
		g g_`i'_`j'_`k'_`y' = g_`i'_`j'_`k' * `y'
		}
		}
		}
		}

		
*** save full data

	* allow fixed effects using id (individual id number)
	xtset id
							
	tempfile aaa
	save `aaa', replace


*===============================================================================
*	Tables 2-4, Figures 3-4
*===============================================================================	

	use `aaa', clear
	
*** Table 2 sample demographics 
	keep if dinner=="CH"   /* can choose any dinner from the panel */
	tab gender
	tab education 
	tab income
	tab age
	
*** Table 3 direction of change in health rating
	use `aaa', clear
	sort dinner condition
	collapse dhr_neg dhr_zero dhr_pos, by(dinner treat)
	outsheet using Table_3.csv, c replace
	
*** Table 4 output: Kolmogorov-Smirnov Tests of equality of the response distributions 
	use `aaa', clear
	set more off 
	foreach i in cal fat sod sug {
		reg i`i' hl tl if dinner=="CH", r
		ksmirnov i`i' if (hl==1 | no==1) & dinner=="CH", by( hl )
		ksmirnov i`i' if (tl==1 | no==1) & dinner=="CH", by( tl )
		ksmirnov i`i' if (hl==1 | tl==1) & dinner=="CH", by( tl )
		} 

*** Figures 3-4 changes in health ratings by dinner, treatment
	use `aaa', clear	
	collapse (mean) dhr adhr (semean) dhr_se=dhr adhr_se=adhr, by(dinner condition)
	g dhr_ub = dhr + dhr_se
	g dhr_lb = dhr - dhr_se
	g adhr_ub = adhr + adhr_se
	g adhr_lb = adhr - adhr_se	
	g treat = 0
	replace treat = 1 if condition=="HL"
	replace treat = 2 if condition=="TL"
	sort dinner treat 
	outsheet dinner condition dhr dhr_lb dhr_ub adhr adhr_lb adhr_ub using Figures_3-4.csv, c replace


*===============================================================================
*    Regressions: Tables 5-6, Figures 6-7
*===============================================================================	
	
*** Regressions all use robust se, for pooled dinners cluster errors at individual id level

	use `aaa', clear

	* set locals for specifications, sugar excluded for relative importance 
	* locals for variable sets	
	local xvars1 "female age45plus colldegree"
	local xvars1a "female x25to34 x35to44 age45to54 age55plus"
	local xdinner "lasag meat stir"
	local xvars2 `xvars1' `xdinner'
	local xvars3 `xvars2' "hsgrad somecoll assocbach graddegree inc25to50 inc50to75 incunder25"
	local imp_rel "ical_rel ifat_rel isod_rel"
	local imptreat_rel "ical_rel_hl ical_rel_tl ifat_rel_hl ifat_rel_tl isod_rel_hl isod_rel_tl"
	local imphrc_rel "ical_rel_hrc ical_rel_hl_hrc ical_rel_tl_hrc ifat_rel_hrc ifat_rel_hl_hrc ifat_rel_tl_hrc isod_rel_hrc isod_rel_hl_hrc isod_rel_tl_hrc"
	local iifemale_rel "ical_rel_female ifat_rel_female isod_rel_female"
	local iiage45_rel "ical_rel_age45plus ifat_rel_age45plus isod_rel_age45plus"
	local iicoll_rel "ical_rel_colldegree ifat_rel_colldegree isod_rel_colldegree"
	local imp_abs "ical_abs ifat_abs isod_abs isug_abs"
	local imptreat_abs "ical_abs_hl ical_abs_tl ifat_abs_hl ifat_abs_tl isod_abs_hl isod_abs_tl isug_abs_hl isug_abs_tl"
	local imphrc_abs "ical_abs_hrc ical_abs_hl_hrc ical_abs_tl_hrc ifat_abs_hrc ifat_abs_hl_hrc ifat_abs_tl_hrc isod_abs_hrc isod_abs_hl_hrc isod_abs_tl_hrc isug_abs_hrc isug_abs_hl_hrc isug_abs_tl_hrc"
	local iifemale_abs "ical_abs_female ifat_abs_female isod_abs_female isug_abs_female"
	local iiage45_abs "ical_abs_age45plus ifat_abs_age45plus isod_abs_age45plus isug_abs_age45plus"
	local iicoll_abs "ical_abs_colldegree ifat_abs_colldegree isod_abs_colldegree isug_abs_colldegree"
	local xdintreat "lasag_hl lasag_tl meat_hl meat_tl stir_hl stir_tl"
	local xvars1treat "female_hl female_tl age45plus_hl age45plus_tl colldegree_hl colldegree_tl"
	local xhrc "hrc hl_hrc tl_hrc"
	local ihrc "hl tl hrc hl_hrc tl_hrc"
	local xvars1hrc "female_hrc female_hl_hrc female_tl_hrc age45plus_hrc age45plus_hl_hrc age45plus_tl_hrc colldegree_hrc colldegree_hl_hrc colldegree_tl_hrc"
	local gamma_abs_num "g_abs_num_cal g_abs_num_fat g_abs_num_sod g_abs_num_sug"
	local gamma_rel_num "g_rel_num_cal g_rel_num_fat g_rel_num_sod"
	local gamma_abs_num_t "g_abs_num_cal_hl g_abs_num_cal_tl g_abs_num_fat_hl g_abs_num_fat_tl g_abs_num_sod_hl g_abs_num_sod_tl g_abs_num_sug_hl g_abs_num_sug_tl"
	local gamma_rel_num_t "g_rel_num_cal_hl g_rel_num_cal_tl g_rel_num_fat_hl g_rel_num_fat_tl g_rel_num_sod_hl g_rel_num_sod_tl"
	local gamma_abs_num_agg "g_abs_num g_abs_num_hl g_abs_num_tl"
	local gamma_rel_num_agg "g_rel_num g_rel_num_hl g_rel_num_tl"
	local dinnertreat "lasag_hl lasag_tl meat_hl meat_tl stir_hl stir_tl"


*** Means and correlations

	tabstat hr dhr adhr hl tl hrc g_abs_num g_rel_num `xdinner' `xvars1' /// 
		, statistics( mean sd p25 p50 p75 ) columns(statistics)
	
	corr hr dhr adhr hrc g_abs_num g_rel_num chick `xdinner' `xvars1'
	
*** Table 5 and Figures 6-7

	foreach z in dhr { 
	foreach i in abs { 	
	
	* model with hrc and gammas, interacted with treatments 
			
		* (1) hrc, treatment dummies, and interactions with treatment 
		reg `z' `ihrc' `gamma_`i'_num_agg' , vce(cluster id)
		est store pool1_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
						
		* (2) same, with individual fixed effects: 
		* note don't really need fe to identify treatment effects given random treatment	
		xtreg `z' `xhrc' `gamma_`i'_num_agg', fe vce(cluster id)
		est store pool1_fe_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
		
		* (3) add dinner dummies to individual fixed effects: 
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner', fe vce(cluster id)
		est store pool1d_fe_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
		
		* (4) include dinner dummies 
		reg `z' `ihrc' `xdinner' `gamma_`i'_num_agg' `xvars1' `xvars1treat', vce(cluster id)
		est store pool3_`i'
						
	* table of results
	
		outreg2 [pool1_`i' pool1_fe_`i' pool1d_fe_`i' pool3_`i'] /// 
			using Table_5.xls, /// 
			sortvar(`ihrc' `gamma_`i'_num_agg' `xdinner' `dinnertreat' `xvars1' `xvars1treat') ///  
			adjr2 excel replace 
		erase Table_5.txt
		
	* the coefficients and ci for all dinners, fe spec 
	* results ordered for plotting interactions
	
		outreg2 [pool1d_fe_`i'] /// 
			using Figures_6-7.xls, /// 
			sortvar(`ihrc' `gamma_`i'_num_agg' `xdinner')  /// 
			adjr2 excel replace stats(ci_high ci_low coef) side noparen noaster
		erase Figures_6-7.txt	
			
		}
		}


*** Table 6	ordered logit 

	use `aaa', clear
	
	*** Ordered logit really only makes sense using hr as dependent var

	* for now we only use numerical nutrient content, which makes the most 
	* sense comparing NO with others
	* model with hrc and gammas, interacted with treatments 
			
	* hrc, treatment dummies, and interactions with treatment
	ologit hr `ihrc' `gamma_rel_num_agg' , vce(cluster id)
	est store ologit1_rel
	ologit hr `ihrc' `gamma_abs_num_agg' , vce(cluster id)
	est store ologit1_abs
												
	* table of results   Table 6
	
	outreg2 [ologit1_rel ologit1_abs] /// 
		using Table_6.xls, /// 
		sortvar(`ihrc' `gamma_`i'_num_agg' `xdinner' `dinnertreat' `xvars1' `xvars1treat' `gamma_`i'_num' `gamma_`i'_num_t') ///  
		excel replace 
	erase Table_6.txt
	

*===============================================================================
*	Additional results, specifications 
*===============================================================================	

*** Balance test of means of demographics and nutritient concern across treatments

	set more off	
	use `aaa', clear
		
	* locals for means tests	
	local vars_cont ical ifat isod isug 	
	local vars_binary female nonwhite age45plus colldegree 
	local vars_all `vars_binary' `vars_cont'
	
	keep `vars_all' id treat
	* keep only one observation per person
	duplicates drop id, force

	* Obtain stats for 2-sided test of means between condition = NO, HL, TL
	* s = , .05, .01 for not sig, 5%, 1%
	
	* means tests
	foreach x in `vars_cont' { 
		quietly sum `x' if treat~=.
		g `x'_m = r(mean)
		g `x'_sd = r(sd)
		quietly ttest `x' if treat==0 | treat==1, by(treat) unequal	
		g `x'_s01 = .
		replace `x'_s01 = .05 if r(p)<=.05
		replace `x'_s01 = .10 if r(p)<=.10			
		g `x'_m0 = r(mu_1)
		g `x'_m1 = r(mu_2)
		quietly ttest `x' if treat==0 | treat==2, by(treat) unequal	
		g `x'_s02 = .
		replace `x'_s02 = .05 if r(p)<=.05
		replace `x'_s02 = .10 if r(p)<=.10			
		g `x'_m2 = r(mu_2)
		quietly ttest `x' if treat==1 | treat==2, by(treat) unequal	
		g `x'_s12 = .
		replace `x'_s12 = .05 if r(p)<=.05
		replace `x'_s12 = .10 if r(p)<=.10	
		}

	* proportions tests
	foreach x in `vars_binary' { 
		quietly sum `x' if treat~=.
		g `x'_m = r(mean)
		g `x'_sd = r(sd)
		quietly prtest `x' if treat==0 | treat==1, by(treat) 	
		g `x'_s01 = .
		replace `x'_s01 = .05 if r(z)>=1.96
		replace `x'_s01 = .10 if r(z)>=1.65
		g `x'_m0 = r(P_1)
		g `x'_m1 = r(P_2)
		quietly prtest `x' if treat==0 | treat==2, by(treat) 	
		g `x'_s02 = .
		replace `x'_s02 = .05 if r(z)>=1.96
		replace `x'_s01 = .10 if r(z)>=1.65
		g `x'_m2 = r(P_2)
		quietly prtest `x' if treat==1 | treat==2, by(treat) 	
		g `x'_s12 = .
		replace `x'_s12 = .05 if r(z)>=1.96
		replace `x'_s01 = .10 if r(z)>=1.65
		}
	
*** create the table
	
	* keep just the first row
	
	keep in 1/1	
	tempfile aaa1
	save `aaa1', replace
	
	* initialize using the overall means and save, then append all the stats
	set more off
	use `aaa1', clear 
	g str stat = "mean" 	
	foreach x in `vars_cont' `vars_binary' { 
		replace `x' = `x'_m		
		}
	keep stat `vars_cont' `vars_binary' 
	tempfile bbb1
	save `bbb1', replace
	
	* loop through other stats	
		
	foreach y in sd m0 m1 m2 s01 s02 s12 { 
		use `aaa1', clear
		g str stat = "`y'" 
	foreach x in `vars_cont' `vars_binary' { 
		replace `x' = `x'_`y'
		} 
		keep stat `vars_cont' `vars_binary' 
		append using `bbb1'
		save `bbb1', replace
		}
	
	* order the rows as desired for table columns		
	g row = 1 if stat=="mean"
	replace row = 2 if stat=="sd"
	replace row = 3 if stat=="m0" 
	replace row = 4 if stat=="m1"
	replace row = 5 if stat=="m2" 
	replace row = 6 if stat=="s01" 
	replace row = 7 if stat=="s02" 
	replace row = 8 if stat=="s12"	
	sort row

	* make a matrix of these reuslts and transpose to get variables in rows
	mkmat row `vars_all', matrix(X)
	matrix Y = X'
	
	* xsvmat allows us to save the row names as a variable
	xsvmat Y, rownames(varname) norestore

	* asterisks for tests of group differences
	foreach x in Y6 Y7 Y8 {
		g str star_`x' = ""
		replace star_`x' = "**" if `x'==float(0.05)
		replace star_`x' = "*" if `x'==float(0.10)
		* float() needed to round binary 
		}
	
	drop if varname=="row"
	keep varname Y1-Y5 star_Y6 star_Y7 star_Y8
		 	
	rename Y1 mean
	rename Y2 sd
	rename Y3 NO
	rename Y4 HL
	rename Y5 TL
	rename star_Y6 HL_NO
	rename star_Y7 TL_NO
	rename star_Y8 TL_HL
	
	outsheet using "$output/balance_test.csv", c replace		
	

*===============================================================================
*	Alternative regression specifications 
*===============================================================================	
	
	set more off	
	use `aaa', clear
	
*** Main results for both relative and absolute measures

	foreach z in dhr { 
	foreach i in rel abs { 	
	
	* model with hrc and gammas, interacted with treatments 
			
		* (1) hrc, treatment dummies, and interactions with treatment 
		reg `z' `ihrc' `gamma_`i'_num_agg' , vce(cluster id)
		est store pool1_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
						
		* (2) same, with individual fixed effects: 
		* note don't really need fe to identify treatment effects given random treatment	
		xtreg `z' `xhrc' `gamma_`i'_num_agg', fe vce(cluster id)
		est store pool1_fe_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
		
		* (3) add dinner dummies to individual fixed effects: 
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner', fe vce(cluster id)
		est store pool1d_fe_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}

		* same as preceding but use log dependent var: 
		xtreg `z'_log `xhrc' `gamma_`i'_num_agg' `xdinner', fe vce(cluster id)
		est store pool1dl_fe_`i'
		
		* add dinner-treatment interactions
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' `dinnertreat', fe vce(cluster id)
		est store pool1dt_fe_`i'
				
		* (4) add demographics to (1)
		reg `z' `ihrc' `gamma_`i'_num_agg' `xvars1' `xvars1treat', vce(cluster id)
		est store pool2_`i'
			* test hl=tl interactions
			test hl_hrc=tl_hrc
			foreach v in g_`i'_num  {
				test `v'_hl = `v'_tl
				}
		
		* include dinner dummies 
		reg `z' `ihrc' `xdinner' `gamma_`i'_num_agg' `xvars1' `xvars1treat', vce(cluster id)
		est store pool3_`i'
						
	* tables of results
	
		outreg2 [pool1_`i' pool1_fe_`i' pool1d_fe_`i' pool1dl_fe_`i' pool1dt_fe_`i' pool2_`i' pool3_`i'] /// 
			using pool_`z'_`i'.xls, /// 
			sortvar(`ihrc' `gamma_`i'_num_agg' `xdinner' `dinnertreat' `xvars1' `xvars1treat') ///  
			adjr2 excel replace 
		erase pool_`z'_`i'.txt
					
		}
		}


*** Run preferred FE specification separately by main demographic groups

	foreach z in dhr { 
	foreach i in rel abs { 	
	
	* dinner dummies with individual fixed effects: 
	
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if female==0, fe vce(cluster id)
		est store male_`i'
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if female==1, fe vce(cluster id)
		est store female_`i'
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if age45plus==0, fe vce(cluster id)
		est store young_`i'
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if age45plus==1, fe vce(cluster id)
		est store old_`i'
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if colldegree==0, fe vce(cluster id)
		est store nodegree_`i'
		xtreg `z' `xhrc' `gamma_`i'_num_agg' `xdinner' if colldegree==1, fe vce(cluster id)
		est store colldegree_`i'
						
	* table of results
	
		outreg2 [male_`i' female_`i' young_`i' old_`i' nodegree_`i' colldegree_`i'] /// 
			using demo_`z'_`i'.xls, /// 
			sortvar(`xhrc' `gamma_`i'_num_agg' `xdinner' ) ///  
			adjr2 excel replace 
		erase demo_`z'_`i'.txt
					
		}
		}
		
*** Separately by dinner
		
	foreach z in dhr {
	foreach i in rel abs { 	
	foreach d in CH LA ML ST { 
				
	* hrc, gamma, demographics  
	
		reg `z' `ihrc' `gamma_`i'_num_agg' `xvars1' `xvars1treat' if dinner=="`d'", vce(r)
		est store adapt_`d'
					
		}

	* tables of results
	
		outreg2 [adapt_CH adapt_LA adapt_ML adapt_ST] /// 
			using adapt_`z'_`i'.xls, /// 
			sortvar(`ihrc' `gamma_`i'_num_agg' `xvars1' `xvars1treat')  /// 
			adjr2 excel replace 
		erase adapt_`z'_`i'.txt

		}
		}
		
*** Absolute change (adhr), pooled dinners

	* For these, including hrc as regressor makes little sense... drop it

	use `aaa', clear 

	foreach z in adhr { 
	foreach i in rel abs { 
	
		* add dinner dummies to individual fixed effects: 
		xtreg `z' `gamma_`i'_num_agg' `xdinner', fe vce(cluster id)
		est store poola_fe_`i'
						
		* demo vars, no fe 
		reg `z' hl tl `xdinner' `gamma_`i'_num_agg' `xvars1' `xvars1treat', vce(cluster id)
		est store poola_`i'
						
	* tables of results
	
		outreg2 [poola_fe_`i' poola_`i'] /// 
			using pool_`z'_`i'.xls, /// 
			sortvar(hl tl `gamma_`i'_num_agg' `xdinner' `xvars1' `xvars1treat' ) ///  
			adjr2 excel replace 
		erase pool_`z'_`i'.txt
		
	* the coefficients and ci for all dinners, spec 2
	* results ordered for plotting interactions
	
		outreg2 [poola_fe_`i'] /// 
			using poola_fe_`z'_`i'_ci.xls, /// 
			sortvar(`gamma_`i'_num_agg' `xdinner')  /// 
			adjr2 excel replace stats(ci_high ci_low coef) side noparen noaster
		erase poola_fe_`z'_`i'_ci.txt	
			
		}
		}	
	
